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noise, and Brownian motion 



Abstract. - We investigate the slow time scales that arise from aging of the paths during the 
process of path aggregation. This is studied using Monte-Carlo simulations of a model aiming 
to describe the formation of fascicles of axons mediated by contact axon-axon interactions. The 
growing axons are represented as interacting directed random walks in two spatial dimensions. 
To mimic axonal turnover, random walkers are injected and whole paths of individual walkers are 
removed at specified rates. We identify several distinct time scales that emerge from the system 
dynamics and can exceed the average axonal lifetime by orders of magnitude. In the dynamical 
steady state, the position-dependent distribution of fascicle sizes obeys a scaling law. We discuss 
our findings in terms of an analytically tractable, effective model of fascicle dynamics. 
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Introduction. — The process of path aggregation is 
a ubiquitous phenomenon in nature. Some examples of 
such phenomena are river basin formation [1] , aggregation 
of trails of liquid droplets moving down a window pane, 
formation of insect pheromone trails [2-4], and of pedes- 
trian trail systems [5,6]. 

Path aggregation has been mathematically studied 
mainly in two classes of models. One of them is known as 
the active- walker models [6] in which each walker in course 
of its passage through the system changes the surround- 
ing environment locally, which in turn influences the later 
walkers. An example of such a process is the ant trail for- 
mation [2,3]. While walking, an ant leaves a chemical trail 
of pheromones which other ants can sense and follow. The 
mechanism of human and animal trail formations is me- 
diated by the deformation of vegetation that generates an 
interaction between earlier and later walkers [6] . A math- 
ematical formalism to study the formation of such trails 
has been developed in Ref. [5,6]. The other class of models 
showing path aggregation deals with non-interacting ran- 
dom walkers moving through a fluctuating environment. 
In Ref. [7] , condensation of trails of particles moving in an 
environment with Gaussian spatial and temporal correla- 
tion is demonstrated analytically. Another example of this 
model class is the Scheidegger river model [1] (and related 
models [8]) which describes the formation of a stream net- 



work by aggregation of streams flowing downhill on a slope 
with local random elevations. 

In this Letter we analyze the dynamics of path aggre- 
gation using a simple model that belongs to the class of 
active walker systems discussed above. The model is simi- 
lar to the one used to study path localization in Ref. [9] . In 
our model, however, we take into account the aging of the 
paths, an important aspect of the active walker models. 
For instance, in ant trail systems, the pheromone trails 
age due to evaporation. In the mammalian trail forma- 
tion the deformations of the vegetation due to the move- 
ment of a mammal decays continuously with time [6]. In 
our model, the individual paths do not age gradually, but 
rather maintain their full identity until they are abruptly 
removed from the system. This particular rule for path 
aging is chosen to allow application of our model to the 
process of axon fasciculation, which we discuss next. 

During the development of an organism, neurons lo- 
cated at peripheral tissues (e.g. the retina or the olfactory 
epithelium) establish connections to the brain via grow- 
ing axons. The growth cone structure at the tip of the 
axon interacts with other axons or external chemical sig- 
nals and can bias the direction of growth when spatially 
distributed chemical signals are present [10]. In the ab- 
sence of directional signals, the growth cone maintains an 
approximately constant average growth direction, while 
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exploring stochastically the environment in the transverse 
direction [11]. The interaction of growth cones with the 
shafts of other axons commonly leads to fasciculation of 
axon shafts [12]. During development a significant portion 
of fully grown neurons die and get replaced by newborn 
neurons with newly growing axons. For certain types of 
neurons (such as the sensory neurons of the mammalian ol- 
factory system) the turnover persists throughout the lifes- 
pan of an animal. In mice, the average lifetime of an olfac- 
tory sensory neuron is 1-2 months [13], which is less than 
one tenth of the mouse lifespan. The mature connectiv- 
ity pattern is fully established only after several turnover 
periods [14]. 

The model we propose in this Letter captures the basic 
ingredients of the process of axon fasciculation, i.e. attrac- 
tive interaction of growth cones with axon shafts, as well 
as neuronal turnover. The main contribution of this Letter 
is a detailed discussion of the slow time scales that emerge 
from the dynamics of our model. Using Monte-Carlo sim- 
ulations we characterize the time scale for the approach 
to steady state and the correlation time within the steady 
state, and show that they can exceed the average axonal 
life time by orders of magnitude. To understand these re- 
sults we formulate an analytically tractable effective single 
fascicle dynamics. This allows us to relate the observed 
slow time scales to the dynamics of the basins of the fasci- 
cles. From the effective fascicle dynamics we derive three 
time scales which we compare to the time scales extracted 
from the Monte-Carlo simulation of the full system. 

For clarity, we stress that the dynamics of our model dif- 
fers substantially from one-dimensional coalescence (A + 
A — > A) [15] or aggregation (mA + nA—> (m + n)A) [16]. 
In our model, there is no direct inter-walker interaction; 
rather, each random walker interacts locally with the trails 
of other walkers. While the stationary properties of the 
system (such as the fascicle size distribution in the steady 
state) may be approximately understood using an analogy 
to one-dimensional diffusion with aggregation, the dynam- 
ical properties (such as the time scale of approach to the 
steady state, and the correlation time in the steady state) 
are undefined in the one dimensional analogy, and require 
understanding based on the full two dimensional model. 
This is in contrast to the situation for path aggregation 
models without turnover: e.g., the Scheidegger river net- 
work model can be mapped onto the Takayasu model of 
diffusion-aggregation in presence of injection of mass in 
one dimension [17]. 

Model and numerical implementation. Each 
growing axon is represented as a directed random walk 
in two spatial dimensions (Fig. la). The random walk- 
ers (representing the growth cones) are initiated at the 
periphery (y = 0, random x) with a birth rate a, and 
move towards the target area (large y) with constant ve- 
locity v y = 1 . In the numerical implementation on a tilted 
square lattice, at each time step the growth cone at (x, y) 
can move to (x — 1, y + 1) (left) or (x + 1, y + 1) (right). 




Fig. 1: (Color online) (a) Interacting directed random walks 
on a tilted square lattice. A random walker (+) represents a 
growth cone. For one walker, the possible future sites (□) and 
their nearest neighbours (o) are marked. The trail of a walker 
(line) models an axon shaft, (b) Typical late-time configuration 
(t = 25T) in a system with L = 800 and No = 100. For the 
fascicle identified at y = 6000 (arrow), D indicates its basin, 
i.e. the interval at the level y = between the right-most 
and left-most axons belonging to the fascicle. The gap E is 
the inter-basin free space at y — 0. Note that the y-coordinate 
cannot be understood as equivalent to time (see the main text). 

(Note that the sites are labeled alternatively by even x 
and odd x at successive y levels.) The probability P{l.r} 
to move left/right is evaluated based on the axon occu- 
pancy at the (x— 1, y + 1) and (x+l,y+l) sites and their 
nearest neighbours (see Fig. la). In the simplest version 
of the model, the interaction is governed by the "always 
attach, never detach" rule: pl = 1 when among the sites 
(x ± 1, y + 1), (x ± 3, y + 1) only (x — 3,y + 1) is already 
occupied; pn = 1 when only (x + 3, y + 1) is occupied; 
Pl = Pr = 1/2 in all other cases. Periodic boundary 
conditions are used in the ^-direction. 

To capture the effect of neuronal turnover, each random 
walker is assigned a lifetime from an exponential distribu- 
tion with mean T. When the lifetime expires, the random 
walker and its entire trail is removed from the system. The 
mean number of axons in the system therefore reaches the 
steady state value N = N exp(—(3y) where -/V = a/ [3, 
and (3 = 1/T is the death rate per axon. In the simula- 
tions, we use T = 10 5 time steps, and restrict our attention 
to y < T/10. The birth rate a is chosen so as to obtain 
the desired number of axons No, or equivalcntly, the de- 
sired axon density p = No/L (p = 1/2 implies an average 
occupancy of one axon per site), where L is the system 
size in x-dircction. The presence of turnover distinguishes 
our model from previous theoretical works on axon fas- 
ciculation [18,19]. The y-coordinate in Fig. la cannot be 
viewed as equivalent to time, and the dynamics at fixed y 
(which is the main focus of this Letter) has no analogy in 
one dimensional models of aggregation or coalescence. 

Mean fascicle size: A typical late-time configuration for 
a system with L = 800 and No = 100 is shown in Fig. lb. 
With increasing fasciculation distance y, the axons aggre- 
gate into a decreasing number m(y) of fascicles. (At a 
given y, two axons are considered to be part of the same 
fascicle if they are not separated by any unoccupied sites.) 
The number of axons in the fascicle is referred to as the 
fascicle size n. The mean fascicle size n at level y may be 
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Fig. 2: (Color online) Approach to the steady state in the 
L = 400, No — 200 system at representative y-levels indicated 
in the legend. The mean fascicle size {n(t;y)} [averaged over 
10 3 initial conditions] approaches n ao {y) as t — * oo. The data 
set labeled as p — 1/8 is from the L = 400, No = 50 system, 
collected at y = 5000. Inset: Power-law growth of mean fascicle 
size rioo—c = 2py b in the steady state. The values of the growth 
exponent and offset are: b = 0.48, c = 1.3 for the L — 800, 
iV = 100 system(x), and b = 0.48, c = 6.9 for the L = 400, 
iV = 200 system(+). 

estimated using the following mean-field argument. Each 
of the m fascicles collects axons that were initiated on 
an interval of length D ~ L/m = Ln/N at the level 
y = (see Fig. 16). The axons initiated at opposite edges 
of the interval are expected to meet within y ~ (D/2) 2 
steps of the random walk in x-direction. Consequently, 
n ~ DN/L ~ 2py 1 / 2 exp(- j3y) for y up to y ~ (L/2) 2 , 
where complete fasciculation (n = N) is expected. Thus 
for y <C (L/2) 2 and /3y <C 1 (which are satisfied in our sim- 
ulations) one obtains the power law growth n ~ 2py 1 ^ 2 . 

Slow time scales. — The measured mean fascicle size, 
obtained by averaging over all the existing fascicles at a 
given y (Fig. 2), grows with time as n — n^ — pexp(— fit) — 
gcxp(— t/r a p), where T ap (y) defines the time scale of ap- 
proach to the steady state value n oa (y). We find that r ap 
is an increasing function of y (up to y ~ (L/2) 2 where a 
single fascicle remains and r ap drops to T) , and can exceed 
the axon lifetime T by orders of magnitude (Figs. 2 and 
36). Asymptotically in y, we find n^ = c + 2py b , with 
6 = 0.480 ± 0.018 (Fig. 2 inset) - in reasonable agreement 
with the mean-field prediction. 

The dynamics in the steady state is characterized by the 
auto-correlation function for the mean fascicle size n(t) at 
a fixed y- level: c(t) = (n(i)n(O)) which fits to the form 
p + qexp(— fit) + rexp(— t/r c ). In Fig. 3a we plot the 
subtracted correlation function g(t) = [c(t) — p]/(q + r). 
The correlation time r c increases with y and significantly 
exceeds the axon lifetime T (Fig. 36). 

Effective fascicle dynamics at fixed y. — We next 
examine the dynamics of individual fascicles. These are 
typically long lived, but at a given y-level the fascicles 
very rarely merge or split (data not shown) . Consequently, 
the number n(t) of axons in each fascicle may be viewed 




Fig. 3: (Color online) (a) Subtracted auto-correlation g(t) at 
y-levels indicated in the legend for the L = 100, No = 50 
system. The time series n(t) is collected over t = 200T to 
2 x 10 4 T; g(t) is further averaged over 30 initial conditions. 
(6) The correlation time r c (□) and approach-to-steady-state 
time scale r ap (o) as a function of y. The theoretical time- 
scales r (filled A), Tf (filled V) (see text) are evaluated from 
the values of a+, b+, c+ in Table I. The solid line is y b with 
b = 1/2. t c (O) for the L = 400, No = 50 system is also shown. 



as given by a stochastic process specified by the rates 
u±(n,y) (for transitions n — > n ± 1). At fixed y, a fas- 
cicle can loose an axon only when the axon dies, thus 
«_ (n) = fin, independent of y. 

The gain rate u+(n,y) is governed by the properties of 
the fascicle basin (see Fig. 16). Under the "always attach, 
never detach" rule, new axons initiated anywhere within 
the basin of size D cannot escape the fascicle. In addition, 
some of the axons born in the two gaps (of size E) flanking 
the basin contribute. Therefore, 

u + = aD/L+(aE/L)[l-U(E,y)] (1) 

where II (£7, y) is the probability that an axon born within 
the gap of size E survives as a single axon at level y. 

The two stochastic variables that fully characterize the 
dynamics of a fascicle are the number of axons n(t) and 
the basin size D{t). In Fig. 4a, we plot n(t) and D(t) 
for a specific fascicle followed over 200T. It is seen that 
n and D tend to co-vary (cross correlation coefficient 
c(D,n) = 0.74). In the following we treat the dynam- 
ics of D as slave to the dynamics of n, i.e. we assume 
that the average separation S = D/{n — 1) between two 
neighbouring axons within the basin is time- independent. 
This implies u + (n) = a + bn with 6 = (3pS. 

The measured average gain rate u+(n) in the steady 
state [20] deviates from linearity at high n, but is well fit 
by u_|_(n) = a + + b + n — c + n 2 (Fig. 46). The quadratic 
correction may be understood as a saturation effect which 
reflects that basins of size D can not exceed 2y or L and 
D > 2y 1 / 2 occur with low probability. The quadratic cor- 
rection to the linear growth of u + with n becomes signifi- 
cant at n > 2y 1 / 2 p sa n. Note that the coefficients a+, b + 
and c+ are functions of y. 

Master equation at fixed y: Let P(n, t) be the distribu- 
tion of fascicles of size n at time t at a fasciculation dis- 
tance y. The master equation of the effective birth-death 
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Fig. 4: (Color online) (a) Time series of number of axons 
n and basin size D for an individual fascicle in a system of 
L — 100 and iVo = 50 at y = 400. The equal time cross 
correlation coefficient averaged over this time-span, c(D, n) = 
[(Dn) — {D)(n)]/ ^J[(D 2 ) - (D) 2 ][(n 2 ) - (n) 2 ] = 0.74. (6) Mean 
gain and loss rates (in units of 0) [averaged over 10 3 initial 
conditions and the interval 100T < t < 150T] as a func- 
tion of fascicle size n. The rates are measured at y = 400 
in a system with L = 100. The fits (lines) are u~ = n, 
= 1.92 + 0.974n - 0.002n 2 (for N = 50), and u {2> = 
1.95 + 0.927n - 0.008n 2 (for iV = 25). 

process may be written as 

P(n,t) = u + {n-l)P(n-l,t) + U-(n+l)P(n + l,t) 
[u+(n)+U-(n)}P(n,t), (2) 

for n > 1. For the boundary state (n = 1) 

P(l,t) = J+(y)+U-{2)P(2,t) - [u + (l)+u-(l)]P(M) 

where J + (y) represents the rate with which new single ax- 
ons appear between existing fascicles at y. In the steady 
state, J+(y) is balanced by the rate with which exist- 
ing fascicles disappear from the system, i.e. J+{y) = 
u_(l)P s (l). The steady state condition P(n,t) = 
then implies pairwise balance, u + {n — l)P s (n — 1) = 
U-(n)P s (n), for all n > 1. Thus the steady state distribu- 
tion is given by P a (n) = J+ (y) Uk=l ^jFj, with the 
normalization condition Yl^Li nP s {n, y) = N. In order 
to obtain a closed- form expression for P s (n, y), we expand 
the pairwise balance condition up to linear order in 1/iV 
and solve to find. 



pP a (n,y) ~ J+{y) n~< exp[-£(n - 1) - n{n - if 



(3) 



where 7 = a+//3 — 1, I = 1 — b+//3 and k = C+/2/3. 

Time scales: Three distinct time scales may be ex- 
tracted from the effective birth-death process. The corre- 
lation time t for the fascicle size n, near the macroscopic 
stationary point n s [u+(n s ) = u_(n s )] can be expressed 
[21] as t = l/(u'_(n s ) - u' + (n s )) = l/(fi - b + + 2c + n s ). 
With a linear approximation of u+(n) = a + + b+n the 
approach-to-steady-state time scale for (n) is r ap = l/(/3 — 
b + ) [21]. We note that the long time scales do not simply 
arise as a consequence of fascicles containing many axons, 
but are due to the dynamics of the fascicle basins. To 
see that, imagine a fascicle with frozen boundary axons, 
for which u+ (a/L)D with D constant. In this case 
u' + (n s ) = and the correlation time r ~ T. To obtain 




Fig. 5: (Color online) Steady state distribution of fascicle sizes 
P s {n,y) [averaged over 10 4 initial conditions and the time in- 
terval 10r < t < 25T] for the iV = 100, L = 800 system 
at y-levels indicated in the legend. Inset: A scaling with 
B — l/(n) and A = in)' 2 ' 1 collapses all data obtained for 
y = 1585, 1995 , 3162 , 5012 , 6310, 7943, 10 4 onto a single curve 
4>{u) = A/"uexp( — vu — Xu 2 ) with u — n/(n) and Af = 274, 
v = 0.78, A = 0.45. 



t > T ', D must co- vary with n. At high y, u' + ps 1/T 
resulting in t ^> T. We note that in the full system, the 
dynamics of the basin size can be viewed as arising from 
the compctction between neighbouring fascicles for basin 
space. 

A third time scale Tf can be defined by J+(y) = m/rf, 
and reflects the rate of turnover of fascicles in the full 
system. Evaluating, Tf = [J ± P s (n,y)dn]/J+(y) using 
P s (n,y) from Eq. 3 (with 7 = 1), we obtain 

T f = (T/2k) \l - (y^e&ie - 2K)erfc(£/2^)) /2 V ^\ . (4) 



Fascicle size distribution in the steady state. 

Following the discussions on effective single fascicle dy- 
namics, we now return to the simulation results for the 
dynamics of the whole system. The steady state is char- 
acterized by the stationary distribution of fascicle sizes 
P s (n,y), defined as the number of fascicles of size n at 
level y. For a system with L = 800 and Nq = 100, P s (n, y) 
is shown at a series of y- levels in Fig. 5. Within the range 
y = 10 3 — 10 4 all data collapse onto a single curve after 
appropriate rescaling (Fig. 5). This data collapse implies 
the scaling law 

P s (n,y) = (n(y))- %1 ^n/(n(y))) (5) 

where the scaling function 4>(u) = Afuexp(—iyu — Xu 2 ). 

Analogy to particle aggregation in one dimen- 
sion. — As we have stressed in the introduction, the full 
dynamics of our system can not be mapped onto the par- 
ticle dynamics of a one-dimensional reaction-diffusion sys- 
tem. At fixed time t (within the steady state), however, 
the aggregation of fascicles with increasing y may be for- 
mally viewed as the evolution of an irreversible aggrega- 
tion process mi + nA — > (m + n)A in one spatial dimen- 
sion, where the y-coordinate (Fig. 1(6)) takes the meaning 



T..A 



Dynamics of path aggregation in the presence of turnover 



y 


a+/(3 


b+/f3 


C+//3 




y 


(S) 


20 


1.96 


0.862 


0.0116 




200 


1.14 


40 


1.94 


0.909 


0.0082 




400 


1.60 


100 


1.96 


0.947 


0.0056 




1000 


2.02 


400 


1.92 


0.974 


0.0022 




10000 


2.21 



Table 1: y-dependence of the fascicle parameters defined in the 
main text. System parameters are L — 100 and No = 50. 

of time. This exhibits analogous scaling properties, but 
with a different scaling function wexp(— Xu 2 ) [16], which 
lacks the exponential part exp(— vu). It is interesting to 
note at this point that to remove the exponential part in 
the expression of the steady state distribution (Eq.3) one 
would require 6+ = j3 which in turn implies r op = oo. In 
other words, without the exponential part exp(— vu) the 
emerging long time scale T ap is lost. This is consistent 
with the fact that the time scales arising from turnover 
are undefined in the one dimensional analogy. 

Scaling of y-dependent parameters. — The expan- 
sion coefficients for u + (n, y), measured at selected y- levels 
in a system with L = 100 and N = 50, are listed in Ta- 
ble I. The scaling property of the distribution of fascicle 
sizes (Eq. 5) implies P s (n,y) = Afy~ 3A b n exp(— vny~ b — 
Xn 2 y~ 2b ). Consistency with Eq. 3 requires J+(y) ~ 
y- 31b , k ~ y~ 2b [i.e. c+ ~ y- 2b ], and t ~ y~ b [i.e. 
(/3 — 6+) ~ y~ b \- This is in reasonable agreement with 
the data in Table I. The coefficient a+ ~ 2(3 for all y re- 
flects 7 ~ 1. Notice that J+(y) oc H(E, y) the probability 
that an axon born in the gap remains single at the level 
y. H(E, y) can be approximated as the survival probabil- 
ity of a random walker moving in between two absorbing 
boundaries (basin borders) each of which is undergoing 
random walk; thus II ~ y~ 3 / 2 [16]. This is consistent with 
the measured exponent in J+(y) ~ y~ 31 °. 

The gradual approach of 6 + w (3pS to (3 with increas- 
ing y may be understood as follows. At low y, fascicles 
are formed preferentially by axons that started with small 
separation S at the y = level, and the typical gap size 
E exceeds S. With increasing y, the smallest gaps are 
removed to become part of fascicles; consequently both E 
and S grow with y, and S approaches 1/p at high y. The 
time averaged S for selected fascicles is shown in the last 
two columns of Table I. 

y -dependence of time scales: The y-dependence of pa- 
rameters discussed above implies the following scaling for 
the theoretically derived time scales: the approach-to- 
steady-state time T ap = l/(/3 — 6+) ~ y b , the correla- 
tion time t ~ T ap ~ y b and the life time of fascicles 
Tf ~ 1/k ~ y 2b . As seen in Fig. 36, the measured 
correlation time r c falls in between the computed time 
scales t and r/, and the three time scales show distinct 
y-dependence. In a system with lower axon density, the 
time scales are reduced (Fig. 36) as the increased inhomo- 
geneity of axon separations at y — leads to a stronger 
deviation of pS from 1. 



Conclusion. — To summarize, we have proposed a 
simple model for axon fasciculation that shows rich dy- 
namical properties. We identified multiple time scales that 
grow with the fasciculation distance y and become 3> T, 
the average lifetime of individual axons. The slow time 
scales do not simply arise as a consequence of fascicles 
containing many axons, but are due to the dynamics of 
the fascicle basins. Our theoretical results have wider rele- 
vance for existing related models (e.g. of insect phcromono 
trails [2-4] and pedestrian trail formation [5,6]), in which 
the slow maturation and turnover of the connectivity pat- 
tern have not been analyzed in detail. 

To conclude, we comment on the applicability of our 
model to the biological system which we described in the 
introduction, i.e. the developing mammalian olfactory sys- 
tem. The model presented in this Letter is readily gener- 
alized to include multiple axon types, type-specific inter- 
actions, and detachment of axons from fascicles [22]. This 
is important as in the olfactory system the nasal epithe- 
lium contains neurons belonging to multiple types, dis- 
tinguished by the expressed olfactory-receptor [12], with 
type specific interaction strengths [23]. (We note, how- 
ever, that preparing transgenic mice expressing only a sin- 
gle olfactory-receptor has recently become a reality [24]). 
A further modification of the basic model that might be 
required is the introduction of history-dependent axonal 
lifetimes [25]. A comparatively simpler task is to relate 
our model to experimental studies of axon growth in neu- 
ronal cell culture [26,27]. 

* * * 

We gratefully acknowledge extensive discussions with 
Paul Feinstein on olfactory development and on the for- 
mulation of our model. 
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